package org.elkastud.model;

/**
 * NeedlemanWunchMatch algorithm for 3 DNASequences implementation.
 * Złożoność czasowa: m*n*o
 * Koszt pamięciowy: m*n*o
 *
 * Implementacja oparta na wyliczaniu macierzy 3-w z oceną dopasowania
 * pochodzącą z jednego z siedmiu kierunków.
 *
 * @author Jakub Skierbiszewski
 */
public class NeedlemanWunchMatch3
{
	/** Porównywane sekwencje (wejściowe) */
	private DNASequence mA, mB, mC;

	/** Macierz oceny */
	private short mMatchArray[][][];

	/** Macierz kierunków, zawiera informację, z którego kierunku pochodzi aktualna ocena */
	private byte mDirections[][][];

	/** Implementacja macierzy podobieństwa poszczególnych nukleotydów */
	private INucleobaseResemblance mResemblance;

	/**
	 * Konstruuje obiekt implementujący algorytm, inicjalizuje pamięć niezbędą do obliczeń
	 * @param a DNASequence
	 * @param b DNASequence
	 * @param c DNASequence
	 */
	public NeedlemanWunchMatch3(DNASequence a, DNASequence b, DNASequence c)
	{
		mResemblance = new NucleobaseSampleResemblance();
		this.mA = a;
		this.mB = b;
		this.mC = c;
		mMatchArray = new short[a.length()+1][b.length()+1][c.length()+1];
		mDirections = new byte[a.length()][b.length()][c.length()];

		short cumulatedGapPenalty = 0;
		for(int i=1; i<a.length()+1 || i<b.length()+1 || i<c.length()+1; ++i)
		{
			cumulatedGapPenalty += mResemblance.getGapPenalty();
			if(i<a.length()) mMatchArray[i][0][0] = cumulatedGapPenalty;
			if(i<b.length()) mMatchArray[0][i][0] = cumulatedGapPenalty;
			if(i<c.length()) mMatchArray[0][0][i] = cumulatedGapPenalty;
		}
	}

	/**
	 * Metoda implementująca algorytm Needlemana-Wuncha dla trzech sekwencji
	 * @return MatchResult3 rezultat dopasowania
	 */
	public MatchResult3 doCalculate()
	{
		//We got 7 directions in which we can go
		short directionValues[] = new short[7];
		byte currDirection;
		Nucleobase empty = new Nucleobase('-');

		//Build matrix
		for(int x=0; x<mA.length()+0; ++x)
		{
			for(int y=0; y<mB.length()+0; ++y)
			{
				for(int z=0; z<mC.length()+0; ++z)
				{
					//Build directions
					for(short i=0; i<7; ++i)
					{
						//Calculate resemblance
						Nucleobase a=empty,b=empty,c=empty;
						if(getX(i)==-1) a=mA.getAt(x);
						if(getY(i)==-1) b=mB.getAt(y);
						if(getZ(i)==-1) c=mC.getAt(z);

						directionValues[i] = (short)
								(mResemblance.getResemblance(a,b,c) +
								 mMatchArray[x+1+getX(i)][y+1+getY(i)][z+1+getZ(i)] );
					}

					//currDirection = getMax(directionValues);
					currDirection = getMaxIndex(directionValues);
					mDirections[x][y][z] = currDirection;
					mMatchArray[x+1][y+1][z+1] = directionValues[currDirection];
				}
			}
		}

		//Build result
		StringBuilder matchX = new StringBuilder();
		StringBuilder matchY = new StringBuilder();
		StringBuilder matchZ = new StringBuilder();

		for(int x=mA.length(), y=mB.length(), z=mC.length(); x>0 && y>0 && z>0; )
		{
			short direction = mDirections[x-1][y-1][z-1];
			if(getX(direction)==0) matchX.append('-');
			else matchX.append(mA.getAt(x-1).getNucleobaseCode());

			if(getY(direction)==0) matchY.append('-');
			else matchY.append(mB.getAt(y-1).getNucleobaseCode());

			if(getZ(direction)==0) matchZ.append('-');
			else matchZ.append(mC.getAt(z-1).getNucleobaseCode());

			x += getX(direction);
			y += getY(direction);
			z += getZ(direction);
		}

		MatchResult3 result = new MatchResult3(
				matchX.reverse().toString(),
				matchY.reverse().toString(),
				matchZ.reverse().toString(),
				mMatchArray[mA.length()][mB.length()][mC.length()]
		);
		return result;
	}

	/**
	 * Sprawdź, czy w podanym kierunku następuje zmiana na pierwszej sekwenji
	 * @param d kierunek
	 * @return	0	nie nastoąpiła zmiana na pierwszej sekwencji
	 *			1	Nastąpiła zmiana na pierwszej sekwencji
	 */
	private byte getX(short d)
	{
		if(d==0 || d==1 || d==2)
			return 0;
		else
			return -1;
	}

	/**
	 * Sprawdź, czy w podanym kierunku następuje zmiana na drugiej sekwenji
	 * @param d kierunek
	 * @return	0	nie nastoąpiła zmiana na drugiej sekwencji
	 *			1	Nastąpiła zmiana na drugiej sekwencji
	 */
	private byte getY(short d)
	{
		if(d==2 || d==5 || d==6)
			return 0;
		else
			return -1;
	}

	/**
	 * Sprawdź, czy w podanym kierunku następuje zmiana na trzeciej sekwenji
	 * @param d kierunek
	 * @return	0	nie nastoąpiła zmiana na trzeciej sekwencji
	 *			1	Nastąpiła zmiana na trzeciej sekwencji
	 */
	private byte getZ(short d)
	{
		if(d==0 || d==3 || d==6)
			return 0;
		else
			return -1;
	}

	/**
	 * Zwraca indeks elementu o najwyższej wartości
	 * @param elements Wyliczone oceny dla siedmiu kierunków
	 * @return index elementu o najwyższej wartości
	 */
	private byte getMaxIndex(short elements[])
	{
		short maxElem = Short.MIN_VALUE;
		byte maxIndex = 0;
		for(byte i=0; i<elements.length; ++i)
		{
			if(maxElem<elements[i])
			{
				maxElem = elements[i];
				maxIndex = i;
			}
		}

		return maxIndex;
	}

}
